{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"color:#777777;background-color:#ffffff;font-size:12px;text-align:right;\">\n",
    "\tprepared by Abuzer Yakaryilmaz (QuSoft@Riga) | November 07, 2018\n",
    "</div>\n",
    "<table><tr><td><i> I have some macros here. If there is a problem with displaying mathematical formulas, please run me to load these macros.</i></td></td></table>\n",
    "$ \\newcommand{\\bra}[1]{\\langle #1|} $\n",
    "$ \\newcommand{\\ket}[1]{|#1\\rangle} $\n",
    "$ \\newcommand{\\braket}[2]{\\langle #1|#2\\rangle} $\n",
    "$ \\newcommand{\\inner}[2]{\\langle #1,#2\\rangle} $\n",
    "$ \\newcommand{\\biginner}[2]{\\left\\langle #1,#2\\right\\rangle} $\n",
    "$ \\newcommand{\\mymatrix}[2]{\\left( \\begin{array}{#1} #2\\end{array} \\right)} $\n",
    "$ \\newcommand{\\myvector}[1]{\\mymatrix{c}{#1}} $\n",
    "$ \\newcommand{\\myrvector}[1]{\\mymatrix{r}{#1}} $\n",
    "$ \\newcommand{\\mypar}[1]{\\left( #1 \\right)} $\n",
    "$ \\newcommand{\\mybigpar}[1]{ \\Big( #1 \\Big)} $\n",
    "$ \\newcommand{\\sqrttwo}{\\frac{1}{\\sqrt{2}}} $\n",
    "$ \\newcommand{\\dsqrttwo}{\\dfrac{1}{\\sqrt{2}}} $\n",
    "$ \\newcommand{\\onehalf}{\\frac{1}{2}} $\n",
    "$ \\newcommand{\\donehalf}{\\dfrac{1}{2}} $\n",
    "$ \\newcommand{\\hadamard}{ \\mymatrix{rr}{ \\sqrttwo & \\sqrttwo \\\\ \\sqrttwo & -\\sqrttwo }} $\n",
    "$ \\newcommand{\\vzero}{\\myvector{1\\\\0}} $\n",
    "$ \\newcommand{\\vone}{\\myvector{0\\\\1}} $\n",
    "$ \\newcommand{\\vhadamardzero}{\\myvector{ \\sqrttwo \\\\  \\sqrttwo } } $\n",
    "$ \\newcommand{\\vhadamardone}{ \\myrvector{ \\sqrttwo \\\\ -\\sqrttwo } } $\n",
    "$ \\newcommand{\\myarray}[2]{ \\begin{array}{#1}#2\\end{array}} $\n",
    "$ \\newcommand{\\X}{ \\mymatrix{cc}{0 & 1 \\\\ 1 & 0}  } $\n",
    "$ \\newcommand{\\Z}{ \\mymatrix{rr}{1 & 0 \\\\ 0 & -1}  } $\n",
    "$ \\newcommand{\\Htwo}{ \\mymatrix{rrrr}{ \\frac{1}{2} & \\frac{1}{2} & \\frac{1}{2} & \\frac{1}{2} \\\\ \\frac{1}{2} & -\\frac{1}{2} & \\frac{1}{2} & -\\frac{1}{2} \\\\ \\frac{1}{2} & \\frac{1}{2} & -\\frac{1}{2} & -\\frac{1}{2} \\\\ \\frac{1}{2} & -\\frac{1}{2} & -\\frac{1}{2} & \\frac{1}{2} } } $\n",
    "$ \\newcommand{\\CNOT}{ \\mymatrix{cccc}{1 & 0 & 0 & 0 \\\\ 0 & 1 & 0 & 0 \\\\ 0 & 0 & 0 & 1 \\\\ 0 & 0 & 1 & 0} } $\n",
    "$ \\newcommand{\\norm}[1]{ \\left\\lVert #1 \\right\\rVert } $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> <font color=\"blue\"> Solutions for </font>Probabilistic Operators</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"task2\"></a>\n",
    "<h3>Task 2</h3>\n",
    "\n",
    "Randomly construct a $ (3 \\times 3 ) $-dimensional probabilistic operator.\n",
    "\n",
    "Randomly determine the entries of the matrix that represents a probabilistic operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Solution</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's start with a zero matrix\n",
    "A = [\n",
    "    [0,0,0],\n",
    "    [0,0,0],\n",
    "    [0,0,0]\n",
    "]\n",
    "\n",
    "# we will randomly construct 3 columns\n",
    "from random import randrange\n",
    "for j in range(3): # each column is iteratively constructed\n",
    "    total = 100 # we will start with 100 and randomly distribute it into three parties\n",
    "    for i in range(2): # we will determine the first two entries randomly \n",
    "        picked_probability = randrange(total)\n",
    "        A[i][j] = picked_probability # the value of the i-th row in j-th column\n",
    "        total = total - picked_probability # remaining part to distribute in the next iteration\n",
    "    A[2][j] = total # the last row (in the j-th column) takes the remaining part after two iterations\n",
    "\n",
    "# let's print matrix A before the normalization\n",
    "# the entries are between 0 and 100\n",
    "print(\"matrix A before normalization:\")\n",
    "for i in range(3):\n",
    "    print(A[i])\n",
    "\n",
    "print(\"the column summations must be 100\")\n",
    "\n",
    "\n",
    "# let's normalize the entries by dividing every element with 100\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        A[i][j] /=100 # shorter form of A[i][j] = A[i][j] / 100\n",
    "        \n",
    "# let's print matrix A after the normalization\n",
    "print() # print an empty line\n",
    "print(\"matrix A after normalization:\")\n",
    "for i in range(3):\n",
    "    print(A[i]) \n",
    "\n",
    "print(\"the column summations must be 1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"task3\"></a>\n",
    "<h3> Task 3 </h3>\n",
    "\n",
    "Write a function in python for Asja's probabilistic operator $ \\mymatrix{cc}{ 0.6 & 0.3 \\\\ 0.4 & 0.7 } $ such that\n",
    "<ul>\n",
    "    <li> it takes a probabilistic state as the input, and </li>\n",
    "    <li> it outputs the new probabilistic state.</li>\n",
    "</ul>\n",
    "   \n",
    "Then, test your function by applying it twice to the starting state $ \\myvector{1 \\\\ 0} $.\n",
    "\n",
    "Remember that: $ \n",
    "    \\myvector{1 \\\\ 0} \\xrightarrow{\\mbox{Asja's coin-flipping protocol}} \\myvector{0.6 \\\\ 0.4}\n",
    "    \\xrightarrow{\\mbox{Asja's coin-flipping protocol}}  \\myvector{0.48 \\\\ 0.52}.\n",
    "$\n",
    "\n",
    "If your function seems to work, then evolve your system for 3, 6, 9, 12, 24, 48, and 96 steps.\n",
    "\n",
    "Is there any pattern?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Solution</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Asja's probabilistic operator (coin flip protocol)\n",
    "A = [\n",
    "    [0.6,0.3],\n",
    "    [0.4,0.7]\n",
    "]\n",
    "\n",
    "# one-step evolution of Asja's probabilistic operator on a given probabilistic state\n",
    "def asja(prelist):\n",
    "    newlist=[0,0]\n",
    "    for i in range(2): # for each row\n",
    "        for j in range(2): # for each column\n",
    "            newlist[i] = newlist[i] + prelist[j] * A[i][j] # summation of pairwise multiplication \n",
    "    return newlist # return the new state\n",
    "            \n",
    "# initial state            \n",
    "state = [1,0]  \n",
    "\n",
    "# after one step\n",
    "state = asja(state)\n",
    "print(\"after one step, the state is\",state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the new state is correct\n",
    "\n",
    "# let's check one more step\n",
    "state = asja(state)\n",
    "print(\"after two steps, the state is\",state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is also correct\n",
    "#\n",
    "# then, let's evolve the system for more steps\n",
    "for i in [3,6,9,12,24,48,96]:\n",
    "    # start from the initial state\n",
    "    state = [1,0]\n",
    "    for t in range(i): # apply asja t times\n",
    "        state = asja(state)\n",
    "    # print the result\n",
    "    print(state,\"after\",(t+1),\"steps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> The system converges to a fixed probabilistic state </b>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"task4\"></a>\n",
    "<h3> Task 4 </h3>\n",
    "\n",
    "Write a function that takes a probabilistic operator and a probabilistic state, and then returns the new probabilistic state.\n",
    "\n",
    "Your function should work for any dimension.\n",
    "\n",
    "Test your function on $ \\mymatrix{ccc}{ 0.4 & 0.6 & 0 \\\\ 0.2 & 0.1 & 0.7 \\\\ 0.4 & 0.3 & 0.3 } $ and $ \\myvector{0.1 \\\\ 0.3 \\\\ 0.6} $. \n",
    "\n",
    "The new probabilistic state should be $ \\myvector{0.22 \\\\ 0.47 \\\\ 0.31} $.\n",
    "\n",
    "Then, evolve your system for 5, 10, 20, and 40 steps.\n",
    "\n",
    "The system should evolve to a fixed probabilistic state.\n",
    "\n",
    "Change your initial state to  $ \\myvector{1 \\\\ 0 \\\\ 0} $, and see whether the converged state is the same or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Solution</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evolve(Op,state):\n",
    "    newstate=[]\n",
    "    for i in range(len(Op)): # for each row\n",
    "        # we calculate the corresponding entry of the new state\n",
    "        newstate.append(0) # we set this value to 0 for the initialization\n",
    "        for j in range(len(state)): # for each element in state\n",
    "            newstate[i] = newstate[i] + Op[i][j] * state[j] # summation of pairwise multiplications\n",
    "    return newstate # return the new probabilistic state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test the function\n",
    "\n",
    "# operator for the test\n",
    "A = [\n",
    "    [0.4,0.6,0],\n",
    "    [0.2,0.1,0.7],\n",
    "    [0.4,0.3,0.3]\n",
    "]\n",
    "\n",
    "# state for test\n",
    "v = [0.1,0.3,0.6]\n",
    "\n",
    "newstate = evolve(A,v)\n",
    "print(newstate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for step in [5,10,20,40]:\n",
    "    new_state = [0.1,0.3,0.6] # initial state\n",
    "    for i in range(step):\n",
    "        new_state = evolve(A,new_state)\n",
    "    print(new_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> The system converges to a fixed probabilistic state </b>.\n",
    "\n",
    "Moreover, the converged probabilistic state is an equal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change the initial state\n",
    "\n",
    "for step in [5,10,20,40]:\n",
    "    new_state = [1,0,0] # initial state\n",
    "    for i in range(step):\n",
    "        new_state = evolve(A,new_state)\n",
    "    print(new_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> The converged probabilistic state does not change with the initial state. </b>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"task7\"></a>\n",
    "<h3> Task 7 </h3>\n",
    "\n",
    "A challenging task.\n",
    "\n",
    "Freivalds reads 50 random strings of length 40. \n",
    "\n",
    "Find the final probabilistic state for each string.\n",
    "\n",
    "Is there any relation between the numbers of $ a $s and $ b $s, say $ N_a $ and $ N_b $, and the probabilities of the first bit being in zero and one, say $ p_0 $ and $ p_1 $?\n",
    "<ul>\n",
    "    <li> When $ N_a > N_b $, then is $ p_0 < p_1 $ or $ p_0 > p_1 $? </li>\n",
    "    <li> When $ N_a < N_b $, then is $ p_0 < p_1 $ or $ p_0 > p_1 $? </li>\n",
    "</ul>\n",
    "\n",
    "Or simply check the signs of $ (N_a - N_b) $ and $ (p_0-p_1) $ for each string.\n",
    "\n",
    "<i>Note that the multiplication of two numbers with the same signs is a positive number, and the multiplication of two numbers with the opposite signs gives a negative number.</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Solution</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for random number generation\n",
    "from random import randrange\n",
    "\n",
    "# we will use evolve function\n",
    "def evolve(Op,state):\n",
    "    newstate=[]\n",
    "    for i in range(len(Op)): # for each row\n",
    "        newstate.append(0)\n",
    "        for j in range(len(state)): # for each element in state\n",
    "            newstate[i] = newstate[i] + Op[i][j] * state[j] # summation of pairwise multiplications\n",
    "    return newstate # return the new probabilistic state\n",
    "\n",
    "# the initial state\n",
    "state = [0.5, 0, 0.5, 0]\n",
    "\n",
    "# probabilistic operator for symbol a\n",
    "A = [\n",
    "    [0.5,  0, 0, 0],\n",
    "    [0.25, 1, 0, 0],\n",
    "    [0,    0, 1, 0],\n",
    "    [0.25, 0, 0, 1]\n",
    "]\n",
    "\n",
    "# probabilistic operator for symbol b\n",
    "B = [\n",
    "    [1, 0, 0,    0],\n",
    "    [0, 1, 0.25, 0],\n",
    "    [0, 0, 0.5,  0],\n",
    "    [0, 0, 0.25, 1]\n",
    "]\n",
    "\n",
    "#\n",
    "# your solution is here\n",
    "#\n",
    "\n",
    "length = 40\n",
    "total = 50\n",
    "# total = 1000 # we will also test our code for 1000 strings\n",
    "\n",
    "# we will check 5 cases\n",
    "# let's use a list \n",
    "cases = [0,0,0,0,0]\n",
    "\n",
    "for i in range(total): # total number of strings\n",
    "    Na = 0\n",
    "    Nb = 0\n",
    "    string = \"\"\n",
    "    state = [0.5, 0, 0.5, 0]\n",
    "    for j in range(length): # generate random string\n",
    "        if randrange(2) == 0: \n",
    "            Na = Na + 1 # new symbol is a\n",
    "            string = string + \"a\"\n",
    "            state = evolve(A,state) # update the probabilistic state by A\n",
    "        else:\n",
    "            Nb = Nb + 1 # new symbol is b\n",
    "            string = string + \"b\"\n",
    "            state = evolve(B,state) # update the probabilistic state by B\n",
    "    # now we have the final state\n",
    "    p0 = state[0] + state [1] # the probabilities of being in 00 and 01\n",
    "    p1 = state[2] + state[3] # the probabilities of being in 10 and 11\n",
    "    #print() # print an empty line\n",
    "    print(\"(Na-Nb) is\",Na-Nb,\"and\",\"(p0-p1) is\",p0-p1)\n",
    "    # let's check possible different cases\n",
    "    \n",
    "    # start with the case in which both are nonzero\n",
    "    # then their multiplication is nonzero\n",
    "    # let's check the sign of their multiplication\n",
    "    if (Na-Nb) * (p0-p1) < 0: \n",
    "        print(\"they have opposite sign\")\n",
    "        cases[0] = cases[0] + 1\n",
    "    elif (Na-Nb) * (p0-p1) > 0: \n",
    "        print(\"they have the same sign\")\n",
    "        cases[1] = cases[1] + 1\n",
    "        \n",
    "    #  one of them should be zero\n",
    "    elif (Na-Nb) == 0:\n",
    "        if (p0-p1) == 0: \n",
    "            print(\"both are zero\")\n",
    "            cases[2] = cases[2] + 1\n",
    "        else: \n",
    "            print(\"(Na-Nb) is zero, but (p0-p1) is nonzero\")\n",
    "            cases[3] = cases[3] + 1\n",
    "    elif (p0-p1) == 0: \n",
    "        print(\"(Na-Nb) is nonzero, but (p0-p1) is zero\")\n",
    "        cases[4] = cases[4] + 1\n",
    "        \n",
    "# check the case(s) observed and the case(s) not observed\n",
    "print() # print an empty line\n",
    "print(cases)        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b> Interpretation </b>\n",
    "\n",
    "Five cases about $ (N_a-N_b) $ and $ (p_0-p_1) $:\n",
    "<ul>\n",
    "    <li> $ case[0] $: they have opposite sign </li>\n",
    "    <li> $ case[1] $: they have the same sign </li>\n",
    "    <li> $ case[2] $: both are zero </li>\n",
    "    <li> $ case[3] $: $ (N_a-N_b) $ is zero, but $ (p_0-p_1) $ is nonzero </li>\n",
    "    <li> $ case[4] $: $ (N_a-N_b) $ is nonzero, but $ (p_0-p_1) $ is zero </li>\n",
    "</ul>\n",
    "\n",
    "We observed only two cases: $ case[0] $ and $ case[2] $.\n",
    "\n",
    "(1) If the numbers of $ a $s and $ b $s are the same, then the probabilities of the first bit being in zero and one are the same.\n",
    "\n",
    "$$\n",
    "    N_a = N_b \\longleftrightarrow p_0 = p_1.\n",
    "$$\n",
    "\n",
    "(2) If the numbers of $ a $s and $ b $s are not the same, then we have only $ (N_a - N_b) \\cdot (p_0-p_1) < 0 $.\n",
    "\n",
    "(2.a) If the number of $ a $s is greater than the number of $ b $s, then the probability of observing $ 0 $ in the first bit ($00$ or $11$) is less than the probability of observing 1 in the first bit ($10$ or $11 $).\n",
    "\n",
    "$$\n",
    "    N_a > N_b \\longrightarrow p_0 < p_1.\n",
    "$$\n",
    "\n",
    "(2.b) If the number of $ a $s is less than the number of $ b $s, then the probability of observing $ 0 $ in the first bit ($00$ or $11$) is greater than the probability of observing 1 in the first bit ($10$ or $11 $).\n",
    "\n",
    "$$\n",
    "    N_a < N_b \\longrightarrow p_0 > p_1.\n",
    "$$\n",
    "\n",
    "<hr>\n",
    "\n",
    "If we have more $ a $s, we expect to observe $ 10 $ and $ 11 $ more than $ 00 $ and $ 01 $.\n",
    "\n",
    "If we have more $ b $s, we expect to observe $ 00 $ and $ 01 $ more than $ 10 $ and $ 11 $.\n",
    "\n",
    "If we have equal numbers of $a$s and $b$s, we expect to observe $ 0 $ and $ 1 $ in the first bit equally."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
